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REVISION 3 CHANGES 


This document is a revised version of the previous one dated June 1972, 
and supersedes that document. The following is a summary of the major techni- 
cal changes included in this revision: 

1. The option of including process noise in the 
cross-track components only has been added 
to the previously existing options of including 
it in all components or not including it at all. 

2. The previously used vehicle switch "s 
and number ,r b M of additionally estimated 
parameters have been combined into one vari- 
able ,f k M . 

3. The magnitude of the process noise to be in- 
cluded is now an input to the routine rather 
than the full 3x3 process noise matrix. 

The matrix is constructed internally from the 
input magnitude. 
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FOREWORD 


This document is one of a series of candidates for inclusion in a future re- 
vision of MSC-04217, "Space Shuttle Guidance Navigation and Control Design Equa- 
tions". The enclosed has been prepared under NAS9- 10268, Task No. 15- A, 
rr GN&C Flight Equation Specification Support ", and applies to function 1 of the 
Orbital Navigation Module (ON2) and function 1 of the Co-orbiting Vehicle Naviga- 
tion Module (ON3) as defined in MSC-03690 Rev. C, "Space Shuttle Orbiter Guid- 
ance, Navigation and Control Software Functional Requirements - Vertical Flight 
Operations", dated 31 July 1972. 



APOLLO Space Guidance Analysis Division 
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NOMENCLATURE 


*d (t> 

b 


c 


nom 


d 



f(q) 

G(t> 

— pole 


_1 

I 


J 


2 


k 


q 


Si 


Perturbing acceleration at time t 

Number of additional quantities, such as landmark 
locations or instrument biases, being estimated 

Constant for adjustment of nominal step-size 

Number of columns in the filter weighting W- matrix 

Primary vehicle covariance matrix (6 x 6) 

Target vehicle covariance matrix (6 x 6) 

Special function of q defined in text 

Gravity gradient matrix 

Unit vector of earth's north polar axis expressed in 
reference coordinates 

Unit vector in the direction of the position vector r 
Three-dimensional identity matrix 

Constant describing dominant term of earth's oblateness 

An index on the filter-weighting matrix elements 

Special function of r and 6 defined in text 

Three-dimensional column vector in the 3x3 process 
noise matrix for either the primary vehicle (i = 3, 4, 5) 
or the target vehicle (i = 9+b, 10+b, 11+b) 
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Qp, Q t 


Process noise matrix (3x3) associated with the primary (P) 
or target (T) vehicle 


-0 

r(t) 

r(t) 


r (t) 
— con 


r (t) 
con 


Geocentric position vector at time tg 

Geocentric position vector at time t 

Magnitude of geocentric position vector 

Reference conic position vector at time t 

Magnitude of reference conic position vector 
at time t 


r 


E 


-F 




s cont 


Mean equatorial radius of the earth 
Geocentric position vector at time tp 
Intermediate values of £ 

Switch indicating whether previous extrapolation is 
to be continued without re-rectification 


S X 

pert 


Switch indicating the perturbing accelerations to 
be included 


s 


q 


s w 


Switch controlling whether process noise is to be 
included in the W-matrix extrapolation , and if so 
whether in all components or only cross -track 
components 

Switch controlling whether state or filter- weighting 
matrix integration is being performed (used only 
internally in routine) 
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t Q Initial time point. Also, time of last rectification 

t F Time to which it is desired to extrapolate (£q, Vq) 

and optionally Wq 

Vq Geocentric velocity vector at time tg 

Vp* Geocentric velocity vector at time tp 

v (t) Reference conic velocity vector at time t 

Wq Filter- weighting matrix at time 

Wp, Filter-weighting matrix at time tp 

w, . Elements of the filter-weighting matrix 

x, 1 

w , . Three-dimensional column vectors into which the 

-k, i 

filter -weighting matrix is partitioned 
x Independent variable in Kepler routine 

x 1 Previous value of x 

£ (t) Vector random variable of dimension b representing 

errors in the additionally estimated quantities such 
as landmark locations or instrument biases 

£^(t) Position deviation vector of true position from 

reference conic position at time t 

6* Magnitude of position deviation vector (temporary 

variable used for rectification test) 

6 Maximum value of I 6 | permitted (used as rectifi- 

max — 

cation criterion) 
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At 


Time-step in numerical integration of differential 
equation 


At 


max 


At 


nom 


€(t) 


Maximum permissible time- step size 

Nominal integration time -step size 

Time convergence tolerance criterion 

Random variable representing error in estimate of 
position vector at time t 


£(t) 


Random variable representing error in estimate of 
velocity vector at time t 


M 

*(t) 


Earth's gravitational parameter 

Velocity deviation vector of true velocity from refer 
ence conic velocity at time t 


Magnitude of velocity deviation vector (temporary 
variable used for rectification test) 


max 


Maximum value of | | permitted (used as rectifi- 

cation criterion) 


Time interval since last rectification 


Previous value of r 


Geocentric latitude 
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1 . 


INTRODUCTION 


The Precision State and Filter Weighting Matrix Extrapola- 
tion Routine provides the capability to extrapolate any spacecraft 
geocentric state vector either backwards or forwards in time through 
a force field consisting of the earth's primary central-force gravita- 
tional attraction and a superimposed perturbing acceleration. The 
perturbing acceleration may be either the single dominant term (J 2 ) 
of the earth's oblateness or a more complete expression involving 
all significant perturbation effects* The Routine also provides the 
capability of extrapolating the filter-weighting matrix along the preci- 
sion trajectory. This matrix, also known as the "W-matrix", is a 
square root form of the error covariance matrix and contains statisti- 
cal information relative to the accuracies of the state vectors and 
certain other optionally estimated quantities. 

On any one call, the routine extrapolates only one state vec- 
tor and only those six rows of the filter- weighting matrix relating to 
this state vector. Two calls are required to extrapolate two separate 
state vectors and a complete filter -weighting matrix pertaining to two 
state vectors. The complete extrapolated filter-weighting matrix is 
obtained by properly adjoining the two separately extrapolated sub- 
matrices of six rows each. 

The routine is merely a coded algorithm for the numerical 
solution of modified forms of the basic differential equations which 
are satisfied by the geocentric state vector of the spacecraft's center 
of mass and by the filter-weighting matrix, namely: 


— r(t> + 
at 2 " 


r 3 (t) 


r(t) = a^t) 


and 

^ W(t) = F(t) W<t) +{jQ[ W T (t)] _1 } 

where a d (t) is the vector sum of all the desired perturbing accelera- 
tions, F(t) is a matrix containing the gravity gradient matrix and the 
identity matrix in its off-diagonal sub-blocks, and Q is the process 
noise matrix. A simplified form of the term in braces is included 
only during phases when process noise is to be introduced into the 
navigation filter to improve the long-term navigation accuracy. 
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Because of its high accuracy and its capability of extrapola- 
ing the filter-weighting matrix, this routine serves as the computa- 
tional foundation for precise space navigation. It suffers from a 
relatively slow computation speed in comparison with the Conic State 
Extrapolation Routine. 
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2 . 


FUNCTIONAL FLOW DIAGRAM 


The Precision State and Filter Weighting Matrix Extrapola- 
tion Routine performs its functions by integrating modified forms of 
the basic differential equations at a sequence of points separated by 
intervals known as time -steps, which are not necessarily of the same 
size. The routine automatically determines the size to be taken at 
each step. 

As shown inFigure 1, a state vector and (optionally) the six 
rows of the filter -weighting matrix relating to this state vector are 
updated one step at a time along the precision trajectory until the 
specified overall transfer time interval is exactly attained. (The size 
of the last time -step is adjusted as necessary to make this possible). 
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3. 


INPUT AND OUTPUT VARIABLES 


The Precision State and Filter Weighting Matrix Extrapolation Routine has 

the following input and output variables: 

Input Variables 

Geocentric state vector to be extrapolated. [ If s Qont f 0, 

(r^, Vq) is last rectified geocentric state vector] 

Time associated with (r Q , v Q ) and . [ If ^ con ^ f 

t Q is last rectification time] 

tp Time to which it is desired to extrapolate (Tq, Vq) and 

optionally 

Spert Switch indicating the perturbing accelerations to be in- 

cluded. < s p er t = 1 implies J 2 oblateness term only; 

Spert > * * m Pli es a more complete perturbing accelera- 
tion model (or models).) 

d Number of columns in filter- weighting matrix <d - 0, 

6, 7, ... , where 0 indicates no W-matrix extrapolation) 

k An index indicating which six rows of the filter-weighting matrix 

are to be extrapolated (hence not used if d = 0). Should be set 
to 0 for the primary vehicle and to 6+b for the target vehicle, 
where b is the number of additionally estimated quantities 
(such as landmark locations or instrument biases) being esti- 
mated. 

W Q Filter- weighting matrix to be extrapolated (optional) 

Sq Switch indicating whether process noise is not to be included 

(8^ - 0), is to be included in all components (s = 1), or is 
to be included in cross -track components only (s^ “ 2) 

q mag Magnitude of the process noise to be included in the W-matrix 

being extrapolated 
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^ Switch indicating whether previous extrapolation is to be 

continued (s cQnt = 1) or not (s cQnt - 0) without re- 
rectification 


(6,y) 

^~con' -con^ 
T 

X 1 

r 1 


Position and velocity deviation vectors 

Conic position and velocity vectors 

Time interval since rectification 

Last value of independent variable in 
Kepler Routine 

Last value of dependent variable in 
Kepler Routine 


At end of previous 
extrapolation [used 

onl y if s cont = 1] 


Output Variables 


(E F'— 


Extrapolated geocentric state vector 


t 


Time associated with (r^,, Vp) and 
[ Will equal tp within tolerance of € t ] 


Wp Extrapolated filter- weighting matrix (only the six rows 

indicated by the index k are extrapolated in a single 
call) 




0 

(6, v) 


(r , v ) 
— con —con 


Last rectified position and velocity vectors 
Time of last rectification 
Position and velocity deviation vectors 
Conic position and velocity vectors 
Time interval since rectification 


For use as input if 

s = 1 on a subse- 
cont 

quent extrapolation 


x 1 Last value of independent variable in Kepler Routine 

T 1 Last value of dependent variable in Kepler Routine 
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‘ 5 ' 'Tb ON OF EQUATIONS 

4 V ;).\ , P^ciad-^, State Extrapolation Equations 

y ’>> perturbing acceleration is small compared with the 

„ -direct numerical integration -of the basics differ- 
hh&f otion ' of the spacecraft state v ect oT^isiinefHc ient . 

£>v£t£;.fd'. o £<;<>h v *v-_o.- uue to Encke is utilized in which only, the devia- 

• a reference conic orbit are numerically integrated* 
les along the reference conic are obtained 



* V - 1 V v „ f • . 1 . 

. ^ the position and velocity vectors, arid / define 

&?.: co;rv orbit, because of the perturbing accelerations*- 

ve? *v\j \y vectors r(t) and v(t) will deviate as "Jr 
s.' * itr *' tf* 1 position and velocity vec-ora r 


on 


(t) 


»mj v a^ ■ i ,f 't-) 'V. cynically extrapolated from r\, and y^. 




i ;> * 


Ci > - £con (t) 


iW = y;{t) - v con i'T) 


fos t|ic, .rector lions, it: ckv? be shown that the position deviation 


t<* 1r 


; *sn£ differ e.^ivA^? equation 


• * 

~~ t rt) + , 

..a 1 " 3 

*3* r 


- | Hq) r(t) + 6(t) 


con 


■?y 


x&sftditions 


6(\ 3 r- ‘s> £(t Q ) = o 


where* 
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q = 


U - 2r) * 6 


2 


r 


f(q) = q 


3 + 3q. + q 2 
1+0+ q) 3 / 2 


and { t ) is the total perturbing acceleration. The above 

second order differential equation in the deviation vector $(t) 
is numerically integrated by a method described in a later sub- 
section. 


The term 


con 


£f(q) r(t) + 6{t)J 


must remain small, i. e, of the same order as a^(t), if the method 
is to be efficient. As the deviation vector Mt) grows in magnitude, 
this term will eventually increase in size. When 


l6(t)| > O.Ollr^t)! or \v(t)\ > 0. 01 1 (t )| 

or when 

I i (t) l> 6 max or l H (t > I > ^max’ 

a new osculating conic orbit is established based on the latest preci- 
sion position and velocity vectors r(t) and v(t), the deviations 6 (t) 
and ^(t) are zeroed, and the numerical integration of 6_(t) and y(t) 
continues. The process of establishing a new conic orbit is called 
rectification. 

The total perturbing acceleration a d ( t ) is in general the 
vector sum of all the desired individual perturbing accelerations com- 
prising the total force field, such as those due to the earth’s oblate- 
ness, the gravitational attractions of the sun and moon, and the earth’s 
atmospheric drag. Since many Shuttle applications will require only 
th^ perturbing effect of the dominant term J 2 of the earth’s oblate- 
ness, the use of only this term has been made a standard option in 
the routine diagrammed in Section 5. However, provision has been 
made for handling a completely general perturbing acceleration. The 
form of this perturbing acceleration will depend primarily upon the 
requirements of the Orbit Navigation function. 
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The explicit expression for the earth's J 2 oblateness accel- 
eration alone is: 


id 


2 


ll j 

2 r 

j2 J 2 [r j 

u 


(1-5 sin 2 0)_i r + 2 sin0_i pole 


where 


i is the unit position vector in reference coordinates, 


— r 


i , is the unit vector of the earth's north polar 
— pole 

axis expressed in reference coordinates, 


sin® = 2 r * 1 pole , 


and 


r^ is the mean equatorial radius of the earth, 
E 


4. 2 Fil ter-Weighting (W) Matrix Extrapolation Equations 

The position and velocity vectors which are maintained by the 
spacecraft's computer are only estimates of the actual values of these 
vectors. As part of the navigation technique it is also necessary for 
the computer to maintain statistical information about the position 
and velocity vectors. Furthermore, in particular applications it is 
necessary to include statistical data on various other quantities, such 
as landmark locations during Orbit Navigation and certain instrument 
biases during Co-orbiting Vehicle Navigation. The filter-weighting 
W-matrix is used for all these purposes. 


If c (t) and t) ( t ) are three dimensional vector random 
variables with zero mean which represent the errors in the estimates 
of a spacecraft's position and velocity at time t, then the six-dimen- 
sional state error covariance matrix E(t) at time t is defined by: 


E(t) 


e ( t ) e ( t ) 


n( t) c( t) 


e (t) T 7 ( t ) 


TJ < t) T) ( t) 


where the bar represents the expected value or ensemble average at 
the fixed time t of each element of the matrix over which it appears. 
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If y ( t ) is a b- dimensional vector random variable with zero 
mean which represents the errors in the estimates of the b additionally 
estimated quantities such as landmark locations or instrument biases, 
then a ( 6 + b) - dimensional state and other parameter covariance matrix 
is defined by: 



€<t) <(t) T 

e(t) t?(t) T 

c(t)V(t) T 

E(t) = 

2<t >e(t) T 

n(t) n(t) T 

T7 (t ) r(t) T 


Y(t) £<t) T 

r(t) n(t) T 

y<t) v(t) T 


Further, if the statistical properties of the positions and ve- 
locities of two separate spacecraft are to be maintained, a twelve- 
dimensional state covariance matrix is defined by: 


E(t) - 


where the subscripts P and T refer to the primary and target ve- 
hicles, respectively* 

And finally, if the statistical properties of the b additionally 
estimated quantities are also to be maintained along with the two state 
vectors, a (12 + b) state and other parameter covariance matrix 
is defined by: 


-p 

€ T 
-P 

T 

£p Up 

T 

lp £t 

T 

— P — T 

-P 

T 

^P 

T 

-P Ip 

„ T 

Ip It 

T 

Up n T 

-T 

T 

S-P 

„ T 
-T 2p 

T 

It £t 

— T 

— T 

ip T 

T 

1 T Ip 

„ T 

H'P S'!* 

_ T 


4-4 



T 

Iplp 

T 

lp Up 

T 

ipi 

T 

— P — T 

_ T 

T 

2p!p 

T 

-p -p 

T 

T 

— P — T 

T 

-P -T 

„ T 
X Ip 

T 

1 Up 

T 

1 1 

T 

X i-T 

T 

2 Hj 

T 

-T-P 

T 

!t Up 

T 

, T 

f rj* f p* 

T 

i-T Hj 

T 

2 T lp 

T 

Up 2p 

T 

Hr x 

T 

H T — T 

T 

2t — T 


Rather than use one of the above covariance matrices in the 
navigation procedure, it is more convenient to use a matrix W (t) 
having the same dimension d as the covariance matrix E (t) and 
defined by: 

E (t) = W (t) W (t) T 


The matrix W (t) is called the filter-weighting matrix, and is in a 

sense a square root of the covariance matrix. 

* 

Extrapolation of the W (t) matrix in time may be made 
by direct numerical integration of the differential equation which it sat- 
isfies . In the one-spacecraft case, this is: 


sr w <*> ■ 

o 

G(t) 

1 3 j °(6 x b) 
O [ 

W(t) + 

0/2) 

IP * 

° o j 

O Q j (S X b) 

°(b X 6) 

i °(bxb) 


( 

^(b x 6) i ^(bxb) 



| ^ 



I 


[ W T (t)] 


(where b - 0, 1, 2, . , . is the number of additionally estimated quantities). 
In the two-spacecraft case, the differential equation is: 

O In 


dT w<11 


G p (t) O 


O 


(b x 6) 


O 

O 


o 

o 


°(6xb) | Q 
I 


O 

O 


I 


o, 


"lb x b) | (b x 6) 


°»««| o 


W(t) + 


The term in braces is added only when process noise is to be included. 
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“O 

° 1 

°(6 x b) 

| O 

O 

\ 

( 

O 


1 0 

O 


< (1/2) 

°(b x 6) 
O 

° 1 

^(b x b) 

| °(b x 6) 

j O 

O 

[ W T (t) ] _1 ,> 

\ 

_o 

0 i 

°{6 x b) 

1 

l O 
1 

Qs - 

/ 


where is the 3x3 identity matrix, the O’s are zero matrices of 
the required dimensions, the G(t ) are the 3x3 conic gravity gradient 
matrices 

G(t) = [3r(t)r(t) T - r 2 (t)I 3 j 

associated with the vehicle under consideration or with the primary (P) 
or target (T) vehicle, and the Q are the likewise associated 3x3 pro- 
cess noise matrices. 

Extrapolation of the W matrix may also be made by the following 
technique, which is somewhat simpler to implement in an on-board computer 
since matrix mainpulations are reduced to more tractable vector manipula- 
tions, and matrix inversion is avoided. 

Let the d x d filter- weighting matrix W - [ w .] be partially 

K» 1 

partitioned into three-dimensional column vectors w, . which bear the 

" * 1 

subscripts of their first component: 







— 0, 0 — 0, 1' ’ —0, 5 

! —0. 6 —0, 5+b 

1 

1 

—0, 6+b — 0, 7+b 

— 0,d-l 

-3,0 — 3, 1* -3, 5 

| w 3> g w 3> 5+b 

l 

1 

j 

— 3, 6+b —3, 7+b 

— 3,d-l 

W 6, 0 w 6, 1 ' * w 6, 1 


‘T 



j W 6, 6 W 6, 5+b 

1 

1 

w 6, 6+b w 6, 7+b 

w 6,d-l 

w 7, 0 w 7, 1 ** w 7, 5 

J W 7, 6 w 7, 5+b 

I till 

1 

1 

1 

1 

w 7, 6+b w 7, 7+b 

w 7.d-l 

W 5+b, 0* • • w 5+b, 5 

[ w 5+b, 6‘ • ’ w 5+b, 5+b 

1 

1 

l 

w 5+b, 6+b W 5+b, 7+b' * 

• w 5+b, d-1 

— 6+b,0' — 6+b, 5 

1 — 6+b, 6’ ’ ’ — 6+b, 5+b 

’ 1 ' 

1 

I 

— 6+b, 6+b -6+b, 7+b‘ ' 

' ' — 6+b, d-1 

— 9+b,0* ‘ *— 9+b, 5 

! — 9+b, 6 - * ‘ — 9+b, 5+6 

1 

i 

1 

— 9+b, 6+b — 9+b, 8+b' 

* ‘ — 9+b, d-1 


* 


The term in braces is added only when process noise is to be included. 
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and let the 3x3 process noise matrices be partitioned into three-dimensional 
column vectors: 

Qp = [ H 3 - £4 • S5 J « Qj = [ Sg+b' 2 . 10+b’ -^ll+b ^ ’ 

Furthermore let the inverse of the filter- weighting matrix be approximated by 
the diagonal matrix 1 (t) whose diagonal elements are the reciprocals of 

diagonal elements of the filter- weighting W- matrix: 




l/w. 


0, 0 




'‘/"d-l.d-l 


Then the previous first order differential equations are equivalent to: 


and 


d t 


2 — 0,i = G(t) ^0,i + (1 / 2w i < i ) ^i 


with 

-3,i 


w 


k, i 

1 2 


dt* 

d 2 

7T 3 " 

with 


f(i - 3, 4, 5 only) \ \ i = 0, 1, . . . (d-1) 


= dT —0, i 

= constant for k>6 


£o,i = G p (t) 


-6+b, i G T (t) -^+b,i + 


w 0ji 4 (1/2 w £ i ) Hi j 
)(i = 3, 4, 5 only) 

(1 ' 2w i.i>ai 

i - 9+b, 10+b, 11+b only)!) i = 0, 1, . . . (d-1) 


w 0 . 
— 3 , 1 


w 


dt — 0, i 


w , 


-9+b, i d t -6+b, i 


w , . - constant for 6 < k < 6+b 

k,i 


The terms in braces are added only when process noise is to be included. 
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(d_1) i = 0, 1 (d-1) 


O 

ii 


When written out in full, the above equations are: 

/ So. i = [ 3 [n< t) -w 0 i ( t) ]i:( t) -r- 2 (t ) w 0 i (t) ] 


, 


d t 


r 5 {t) 


+ (1 ^ 2w i,i ) Hi (i = 3,4, 5 only) 


with 


w , . 

— 3, i 


W k,i 

and 

j 2 


dX 


constant for k ^ 6 


d t 


2 -0, i 


~T^-[ 3 [lp(t)* w 0i (t)Jr p (t) - r p 2 (t)w 0 .(t)j 


| (1/2 i )^. (i = 3, 4, 5 only) 


d t 


— 6+b, i 


r x 5 (t) 


[ 3 [, 


x (t) ■ Wg + ^j(t) £ x (t) - r x (t) 


S 6+ b.i (t) ] 


| U / 2w i, i>Si (i = 9+ t>, 10+b. 11+b only) | 


with 


— 3, i 

— 9+b, i 


w 


k, i 


TT ^o,i 

d 

dt —6+b, i 

constant for 6<k < 6 + b 


These second-order differential equations may be integrated using the 
same numerical integration technique as is used for the spacecraft 
position vector. The vectors w 3 . and w g+b . bear the same relation- 
ship to the spacecraft velocity vector as the vectors w n . and w„ 

0 * i 6+b, i 

bear to the spacecraft position vector, and w ~ . ar)d w_ , . area 

— o,\ — 9+b, l 

by-product of the numerical integration of w . and w e iust as 
the velocity vector is a by-product of the numerical integration of the 
position vector. 


The terms in braces are added only when process noise is to be included. 


4-8 



4. 3 Numerical Integration Method 

The extrapolation of inertial state vectors and filter weight- 
ing matrices requires the numerical solution of two second -order 
vector differential equations, which are special cases of the general 
form 

j2 


dt 


where 


£<t) = f(t, Z <t), z(t) ) 


z=—y_ 

dt 


Nystrom’s standard fourth-order method is utilized to numerically 
solve this equation. The algorithm for this method is: 


in+1 = + ^ + k 3 > ( At) 2 

b 

inH-5n + ^<!‘l + 2 ‘‘2 + 2 t3 + !S4> At 


where 


h. = f - (t n * 2n> 


—2 


—3 


l - (, n + in * J Jn At + jV + “ » 

L (t n +At - Z n + Jn A * t jS 3 ( A ‘> 2 . Jn + iSa at > 

Zv * ^ - z(t n > 


and 


Vi • v + ^ 
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As can be seen, the method requires four evaluations of 
f (t, z ) per integration step At as does the classical fourth-order 
Runge-Kutta method when it is extended to second -order equations. 
However, if f is independent of z, then Nystrom's method above only 
requires three evaluations per step since = kg. (Runge-Kutta f s 
method will still require four). 


The integration time step At may be varied from step to 
step. The nominal integration step size is 


At - c r ~/ n/ m 

nom nom con 




where c is a program constant. (The value c = 0. 3 is 

nom nom 

recommended and implies that about 21 steps will be taken per trajec- 
tory revolution). The actual step-size is however limited to a maxi- 
mum of At , which is also a program constant. (A value of about 
max 

4000 seconds is suggested.) Also, in the last step, the actual step 
size is taken to be the interval between the end of the previous step 
and the desired integration endpoint, so that the extrapolated values 
of the state or W-matrix are immediately available. Thus the integra- 
tion step-size At is given by the formula 


At = + minimum (ftp- t|, At^, At^) 

where t^is the desired integration end-point and t is the time at the 
end of the previous step. The plus sign is used if forward extrapola- 
tion is being performed, while the negative sign is used in the back- 
dating case. 
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DETAILED FLOW DIAGRAMS 


This section contains detailed flow diagrams of the Precision State and 
Filter Weighting Matrix Extrapolation Routine. 

Each input and output variable in the routine and subroutine call state- 
ments can be followed by a symbol in brackets. This symbol identifies the 
notation for the corresponding variable in the detailed description and flow dia- 
grams of the called routine. When identical notation is used, the bracketed 
symbol is omitted. 
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UNIVERSAL 

CONSTANTS 



PROGRAM 

CONSTANTS 


e. . c , At , 
t nom max 

^ max J v max 


INPUT VARIABLES 
HO'Xo’W s pert ; 
d, k. W Q ; 


s q ' ^mag ' 

s . , 6 , v , r . v , 1 ”, x'.t’ 
cont — — —con —con 


h_ = (0, 0, 0) 

V = (0, 0, 0) 


r = r n 
—con —0 


Xcon = ^0 


= 0, x’ = 0 

= 0 , T' = 0 


r 

= r +6 

— 0 

—con — 

V n 

= v +v 

-0 

— con — 

in 

- t + T 

0 

0 


6 , Ucpnl or>6 

100 max 


\ 100 ma> 



r 3 / 2 /jT 
nom con / v 

sign ( t F - t ) 

a min fl t - t|. At At l 

|_l r \ nom max J 




r 

c: 

r 

+ <5 

— F 


— con 


>1 


—con 

+ V 

W F 

= 

~k,0 




— k+3. 

0 — ] 


OUTPUT 

VARIABLES f 

ilp'-Zp’ ■ 

— O' — O' t Q '— ■ - '- con’-con' t . x '. t ' 


Figure 2a. Detailed Flow Diagram 











Figure 2b. Detailed Flow Diagram 
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- a + ^ + "g(kj + k2 + k3)AtjAt 


v - + -g (Icj^ -»- 2kg + 2 kg + k ^ 













6. SUPPLEMENTARY INFORMATION 


Encke's technique is a classical method in astrodynamics 
and is described in all standard texts, for example Battin (1964). 

The f (q) function used in Encke's technique (and in the lunar-solar 
perturbing acceleration computations) has generally been evaluated 
by a power series expansion; the closed form expression given here 
was derived by Potter, and is described in Battin (1964). 

The oblateness acceleration in terms of a general spherical 
harmonic expansion may be calculated in a variety of ways; three 
different recursive algorithms are given in Gulick (1970). For low 
order expansions, especially those involving mostly zonal terms, an 
explicit formulation is generally superior computation-time-wise, as 
only the non-zero terms enter into the calculation. The general ex- 
pression for the zonal terms is given by Battin (1964), while Zeldin 
and Robertson (1970) give explicit analytic expressions for each of 
the tesseral terms up through fifth order; hence all combinations of 
terms may easily be included in the oblateness acceleration by con- 
sulting the formulations in these references. 

A full discussion of the use of covariance matrices in space 
navigation is given in Battin (1964). Potter (196 3) suggested the use of 
the W -matrix and developed several of its properties. It should be noted 
that strictly the gravity gradient matrix G(t) should also include the 
gradient of the perturbing acceleration; however, these terms are so 
small that they may be neglected for our purposes. The use of only 
the conic gravity gradient, however, does not imply the W- matrix is 
being extrapolated conically. (Conic extrapolation of the W-matrix 
can be performed by premultiplying the W-matrix by the conic state 
transition matrix, which can be expressed in closed form). Rather 
the W-matrix is here extrapolated along the precision (perturbed ) 
trajectory, as can be seen from the detailed flow diagram of Section 
5. 

The expression for the inclusion of process' noise in the 
differential equation satisfied by the filter -weighting matrix is taken 
from Gustafson and Kriegsman (1970), page 7. 
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The Nystrom numerical integration technique was first con- 
ceived by Nystrom ( 1925 ), and is described in all standard texts on 
the numerical integration of ordinary differential equations, such as 
Henrici (1962), Parametric studies carried out by Robertson (1970) 
on the general fourth -order Runge-Kutta and Nystrom integration 
techniques indicate that the "classic” techniques are the best overall 
techniques for a variety of earth orbiting trajectories in the sense of 
minimizing the terminal position error for all the trajectories, 
although for any one trajectory a special technique can general^ be 
found which decreases the position error after ten steps by one or 
two orders of magnitude for only that trajectory. The classical 
fourth -order Runge-Kutta and Nystrom techniques are approximately 
equally accurate, but the latter possesses the computational advant- 
age of requiring one less perturbing acceleration evaluation per step 
when the perturbing acceleration is independent of the velocity. This 
fact has been taken into account in the detailed flow diagram of Section 
5, in that the extra evaluation is performed only when the perturbing 
acceleration depends explicitly on the velocity. Some past Apollo ex- 
perience has suggested that extra evaluation effect with drag is so 
small as to be negligible; further analysis will confirm or deny this 
for the Space Shuttle. In regard to step-size, the constants and the 
functional form of the nominal and maximum time ^step expressions 
have been determined by Marscher (1965). 
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